################################################################################################
# This script takes the clean recoded SV data and uses the histogram from the Population Imm. 
# proposal to bootstrap the SV exp data. The reasoning behind this bootsrap is because it is  
# plausible that SV immigration proposals are a distinctly different sample
# from that observed in the QV immigration, Bootstrap will sample with replacement. 
# 
# It goes through 10000 bootstraps conditonal on imitating the Pop Imm Pref histogram and runs 
# SV Rules that would be used in the normal bootstrap for QV/SV. It saves the results of the
# rules and also a copy of the distribution preferences for each proposal in each bootstrap 
# iteration. Exports these files locally to the working directory.
# 
# This script takes the clean SV data and bootstraps  the exp data to determine welfare under 
# different rules to cast the BV. The following rules are tested:
# 
# - Rule A: Cast BV on issue subject chose in experiment
# - Rule B: Cast BV according to the strategies and probabilities estimates in the stat model.
# - Rule C: Cast BV on an issue with most points.
# - Rule D: Cast BV on an issue with most points or at random (decided 50/50 chance).
# - Rule E: Cast BV on an issue at random.
# 
# Rewrites the BV columns based on these and then calculates some welfare measures. 
# 
# All data is the saved locally to a dataframe and then exported to the folder
# SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/
# 
# Author: Luis S.
# Last Modified: 7.25.16
################################################################################################

library(dplyr)

setwd("C:/Users/Luis/Dropbox/Research/SV vs QV California Experiment/Experiment Analysis & Results")

SV_data <- read.csv("California Stage 2 - SV & QV Data/Clean Data/California_MTurk_SV_Stage_2_clean_prefrecoded.csv", header = T, stringsAsFactors = F)
QV_data <- read.csv("California Stage 2 - SV & QV Data/Clean Data/California_MTurk_QV_Stage_2_clean_prefrecoded.csv", header = T, stringsAsFactors = F)


###################################
#### SV Vote & Welfare Summary ####
##################################%

#--- Vote Summary and Welfare Summary Fns

SV_VoteSummary <- function(tempBootstrap)
{
  #--- Helper function for SV Summary
  countBV <- function(tempBV, SV_Summary, Proposal, SV)
  {
    if(SV == 1)
    {
      tempAbs <- filter(tempBV, SV1_Vote == Proposal, SV1_VotePosition == "Abstain")$BV1
      tempAgainst <- filter(tempBV, SV1_Vote == Proposal, SV1_VotePosition == "Against")$BV1
      tempInFav <- filter(tempBV, SV1_Vote == Proposal, SV1_VotePosition == "InFavor")$BV1
    }
    
    if(SV == 2)
    {
      tempAbs <- filter(tempBV, SV2_Vote == Proposal, SV2_VotePosition == "Abstain")$BV2
      tempAgainst <- filter(tempBV, SV2_Vote == Proposal, SV2_VotePosition == "Against")$BV2
      tempInFav <- filter(tempBV, SV2_Vote == Proposal, SV2_VotePosition == "InFavor")$BV2
    } 
    
    SV_summary[[paste(Proposal,"_onlyBV",SV,sep="")]][1] <- ifelse(length(tempAbs) > 0, tempAbs, 0)
    SV_summary[[paste(Proposal,"_onlyBV",SV,sep="")]][2] <- ifelse(length(tempAgainst) > 0, tempAgainst, 0)
    SV_summary[[paste(Proposal,"_onlyBV",SV,sep="")]][3] <- ifelse(length(tempInFav) > 0, tempInFav, 0)
    
    return(SV_summary)
  }
  
  
  detEffResult <- function(Position, Intensity)
  {
    # synthesize into a single dataframe and separate by positions
    temp <- data.frame(RefPosition = Position, PrefIntensity = Intensity, stringsAsFactors = F)
    
    # separate by infavor and against
    tempInFav <- temp %>% filter(RefPosition == "InFavor") 
    tempAgainst <- temp %>% filter(RefPosition == "Against")
    
    # compare intensity of preferences and determine winner
    tempInFav <- sum(tempInFav$PrefIntensity)
    tempAgainst <- sum(tempAgainst$PrefIntensity)
    effWinner <- ifelse(tempInFav == tempAgainst, sample(c("InFavor","Against"), 1, prob = c(.5,.5), replace = T),
                        ifelse(tempInFav > tempAgainst, "InFavor", "Against"))
    
    return(effWinner)
  }
  
  
  # #FOR TESTING ONLY
  # tempBootstrap <- BS_data
  # tempBootstrap <- expData
  
  #--- Determine the direction of the bonus vote. 
  tempBootstrap <- mutate(tempBootstrap, SV1_VotePosition = NA, SV2_VotePosition = NA)
  
  for(row in 1:nrow(tempBootstrap))
  {
    # Mutate dataframe with bonus vote position
    if(tempBootstrap$SV1_Vote[row] == "BilingualEduc") {tempBootstrap$SV1_VotePosition[row] <- tempBootstrap$Ref1_BilingualEduc[row]}
    if(tempBootstrap$SV1_Vote[row] == "TeacherTenure") {tempBootstrap$SV1_VotePosition[row] <- tempBootstrap$Ref2_TeacherTenure[row]}
    if(tempBootstrap$SV1_Vote[row] == "PublicBonds") {tempBootstrap$SV1_VotePosition[row] <- tempBootstrap$Ref3_PublicBonds[row]}
    if(tempBootstrap$SV1_Vote[row] == "Immigration") {tempBootstrap$SV1_VotePosition[row] <- tempBootstrap$Ref4_Immigration[row]}
    
    if(tempBootstrap$SV2_Vote[row] == "BilingualEduc") {tempBootstrap$SV2_VotePosition[row] <- tempBootstrap$Ref1_BilingualEduc[row]}
    if(tempBootstrap$SV2_Vote[row] == "TeacherTenure") {tempBootstrap$SV2_VotePosition[row] <- tempBootstrap$Ref2_TeacherTenure[row]}
    if(tempBootstrap$SV2_Vote[row] == "PublicBonds") {tempBootstrap$SV2_VotePosition[row] <- tempBootstrap$Ref3_PublicBonds[row]}
    if(tempBootstrap$SV2_Vote[row] == "Immigration") {tempBootstrap$SV2_VotePosition[row] <- tempBootstrap$Ref4_Immigration[row]}
  }
  
  
  # assign("Testing_BootstrappedSample", tempBootstrap, envir=.GlobalEnv)
  
  # Determine count w/o bonus votes #***#
  tempRef1 <-  c(sum(tempBootstrap$Ref1_BilingualEduc == "Abstain"), sum(tempBootstrap$Ref1_BilingualEduc == "Against"), sum(tempBootstrap$Ref1_BilingualEduc == "InFavor"))
  tempRef2 <-  c(sum(tempBootstrap$Ref2_TeacherTenure == "Abstain"), sum(tempBootstrap$Ref2_TeacherTenure == "Against"), sum(tempBootstrap$Ref2_TeacherTenure == "InFavor"))
  tempRef3 <-  c(sum(tempBootstrap$Ref3_PublicBonds == "Abstain"), sum(tempBootstrap$Ref3_PublicBonds == "Against"), sum(tempBootstrap$Ref3_PublicBonds == "InFavor"))
  tempRef4 <-  c(sum(tempBootstrap$Ref4_Immigration == "Abstain"), sum(tempBootstrap$Ref4_Immigration == "Against"), sum(tempBootstrap$Ref4_Immigration == "InFavor"))
  
  SV_summary <- data.frame(Position = c("Abstain", "Against", "InFavor"), BilingualEduc_NoBV = tempRef1,
                           TeacherTenure_NoBV = tempRef2, PublicBonds_NoBV = tempRef3, Immigration_NoBV = tempRef4)
  
  
  # Determine count w/ SV1 & SV2 bonus votes
  tempBV1 <- group_by(tempBootstrap, SV1_Vote, SV1_VotePosition) %>% summarise(BV1 = n())
  
  SV_summary <- countBV(tempBV1, SV_summary, "BilingualEduc", SV = 1)
  SV_summary <- countBV(tempBV1, SV_summary, "TeacherTenure", SV = 1)
  SV_summary <- countBV(tempBV1, SV_summary, "PublicBonds", SV = 1)
  SV_summary <- countBV(tempBV1, SV_summary, "Immigration", SV = 1)
  
  tempBV2 <- group_by(tempBootstrap, SV2_Vote, SV2_VotePosition) %>% summarise(BV2 = n())
  
  SV_summary <- countBV(tempBV2, SV_summary, "BilingualEduc", SV = 2)
  SV_summary <- countBV(tempBV2, SV_summary, "TeacherTenure", SV = 2)
  SV_summary <- countBV(tempBV2, SV_summary, "PublicBonds", SV = 2)
  SV_summary <- countBV(tempBV2, SV_summary, "Immigration", SV = 2)
  
  # Clean up the SV_Summary dataframe
  SV_summary <- SV_summary %>%
    mutate(BilingualEduc_withBV1 = BilingualEduc_NoBV + BilingualEduc_onlyBV1,
           BilingualEduc_withBV2 = BilingualEduc_NoBV + BilingualEduc_onlyBV2,
           TeacherTenure_withBV1 = TeacherTenure_NoBV + TeacherTenure_onlyBV1,
           TeacherTenure_withBV2 = TeacherTenure_NoBV + TeacherTenure_onlyBV2,
           PublicBonds_withBV1 = PublicBonds_NoBV + PublicBonds_onlyBV1,
           PublicBonds_withBV2 = PublicBonds_NoBV + PublicBonds_onlyBV2,
           Immigration_withBV1 = Immigration_NoBV + Immigration_onlyBV1,
           Immigration_withBV2 = Immigration_NoBV + Immigration_onlyBV2) %>%
    select(-matches("only"))
  
  assign("Testing_SV_Summary_Pre", SV_summary, envir=.GlobalEnv)
  
  SV_summary <- SV_summary[-1,]
  
  row.names(SV_summary) <- seq(1,nrow(SV_summary),1)
  
  SV_summary <- as.data.frame(t(SV_summary), stringsAsFactors = F)
  SV_summary <- tibble::rownames_to_column(SV_summary, "Position")
  colnames(SV_summary) <- as.character(SV_summary[1,])
  SV_summary <- SV_summary[-1,]
  
  # Calculate desired values for each of the referenda
  SV_summary <- arrange(SV_summary, Position)
  
  for(ref in 1:4) # iterates in multiples of four to get data for each referenda
  {
    index_noBV <- (ref*3)-2
    index_BV1 <- (ref*3)-1
    index_BV2 <- (ref*3)
    
    # Calculating simple majority result w/o BV (rows 1,4,7,10) & Determine the simple majority winner
    if(SV_summary$InFavor[index_noBV] == SV_summary$Against[index_noBV])
    {
      SV_summary$NoBV_Majority[index_noBV] <- "Tie"
      SV_summary$NoBV_Majority[index_BV1] <- "Tie"
      SV_summary$NoBV_Majority[index_BV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      SV_summary$NoBV_Result[index_noBV] <- tieBreaker
      SV_summary$NoBV_Result[index_BV1] <- tieBreaker
      SV_summary$NoBV_Result[index_BV2] <- tieBreaker
      
    } else if(SV_summary$InFavor[index_noBV] > SV_summary$Against[index_noBV])
    {
      SV_summary$NoBV_Majority[index_noBV] <- "InFavor"
      SV_summary$NoBV_Majority[index_BV1] <- "InFavor"
      SV_summary$NoBV_Majority[index_BV2] <- "InFavor"
    } else {
      SV_summary$NoBV_Majority[index_noBV] <- "Against"
      SV_summary$NoBV_Majority[index_BV1] <- "Against"
      SV_summary$NoBV_Majority[index_BV2] <- "Against"
    }
    
    SV_summary <- mutate(SV_summary, NoBV_Result = ifelse(NoBV_Majority == "Tie", NoBV_Result, NoBV_Majority))
    
    
    # Calculating BV1 result (rows 2,5,8,11) & Determine the BV1 winner
    if(SV_summary$InFavor[index_BV1] == SV_summary$Against[index_BV1])
    {
      SV_summary$BV1_Majority[index_noBV] <- "Tie"
      SV_summary$BV1_Majority[index_BV1] <- "Tie"
      SV_summary$BV1_Majority[index_BV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      SV_summary$BV1_Result[index_noBV] <- tieBreaker
      SV_summary$BV1_Result[index_BV1] <- tieBreaker
      SV_summary$BV1_Result[index_BV2] <- tieBreaker
      
    } else if(SV_summary$InFavor[index_BV1] > SV_summary$Against[index_BV1])
    {
      SV_summary$BV1_Majority[index_noBV] <- "InFavor"
      SV_summary$BV1_Majority[index_BV1] <- "InFavor"
      SV_summary$BV1_Majority[index_BV2] <- "InFavor"
    } else {
      SV_summary$BV1_Majority[index_noBV] <- "Against"
      SV_summary$BV1_Majority[index_BV1] <- "Against"
      SV_summary$BV1_Majority[index_BV2] <- "Against"
    }
    
    SV_summary <- mutate(SV_summary, BV1_Result = ifelse(BV1_Majority == "Tie", BV1_Result, BV1_Majority))
    
    # Calculating BV2 result (rows 2,5,8,11) & Determine the BV2 winner
    if(SV_summary$InFavor[index_BV2] == SV_summary$Against[index_BV2])
    {
      SV_summary$BV2_Majority[index_noBV] <- "Tie"
      SV_summary$BV2_Majority[index_BV1] <- "Tie"
      SV_summary$BV2_Majority[index_BV2] <- "Tie"
      
      tieBreaker <- ifelse(sample(2, 1, replace=T) > 1, "InFavor", "Against")
      
      SV_summary$BV2_Result[index_noBV] <- tieBreaker
      SV_summary$BV2_Result[index_BV1] <- tieBreaker
      SV_summary$BV2_Result[index_BV2] <- tieBreaker
      
    } else if(SV_summary$InFavor[index_BV2] > SV_summary$Against[index_BV2])
    {
      SV_summary$BV2_Majority[index_noBV] <- "InFavor"
      SV_summary$BV2_Majority[index_BV1] <- "InFavor"
      SV_summary$BV2_Majority[index_BV2] <- "InFavor"
    } else {
      SV_summary$BV2_Majority[index_noBV] <- "Against"
      SV_summary$BV2_Majority[index_BV1] <- "Against"
      SV_summary$BV2_Majority[index_BV2] <- "Against"
    }
    
    SV_summary <- mutate(SV_summary, BV2_Result = ifelse(BV2_Majority == "Tie", BV2_Result, BV2_Majority))
  }
  
  # #FOR TESTING ONLY
  # testingSummary <- SV_summary
  
  SV_summary <- SV_summary[grep("NoBV", SV_summary$Position), ]
  
  SV_summary$Position <- gsub("_NoBV", "", SV_summary$Position)
  
  assign("Testing_SV_Summary_Mid", SV_summary, envir=.GlobalEnv)
  
  # Clean summary table
  
  SV_summary <- select(SV_summary, Position, matches("Result")) %>% rename(Referenda = Position) %>% mutate(Eff_Result = NA)
  
  # Determine the efficient outcome and add it to SV Summary
  
  SV_summary$Eff_Result[SV_summary$Referenda == "BilingualEduc"] <- detEffResult(tempBootstrap$Ref1_BilingualEduc, tempBootstrap$EducationPref)
  SV_summary$Eff_Result[SV_summary$Referenda == "Immigration"] <- detEffResult(tempBootstrap$Ref4_Immigration, tempBootstrap$ImmigrationPref)
  SV_summary$Eff_Result[SV_summary$Referenda == "PublicBonds"] <- detEffResult(tempBootstrap$Ref3_PublicBonds, tempBootstrap$BondsPref)
  SV_summary$Eff_Result[SV_summary$Referenda == "TeacherTenure"] <- detEffResult(tempBootstrap$Ref2_TeacherTenure, tempBootstrap$TeachersPref)
  
  # Determine if the BV helped the minority win
  
  SV_summary <- SV_summary %>%
    mutate(SV1_Winner = ifelse(NoBV_Result == BV1_Result, "Majority", "Minority"),
           SV2_Winner = ifelse(NoBV_Result == BV2_Result, "Majority", "Minority")) %>%
    mutate(IsEff_NSV = ifelse(NoBV_Result == Eff_Result, "Efficient", "Inefficient"),
           IsEff_SV1 = ifelse(BV1_Result == Eff_Result, "Efficient", "Inefficient"),
           IsEff_SV2 = ifelse(BV2_Result == Eff_Result, "Efficient", "Inefficient"))
  
  # assign("Testing_SV_Summary_Post", SV_summary, envir=.GlobalEnv)
  assign("SV_summary", SV_summary, envir=.GlobalEnv) # needed for later... do not remove!!
  
  NSV_Results <- SV_summary$NoBV_Result
  Eff_Results <- SV_summary$Eff_Result
  SV1_Results <- SV_summary$BV1_Result
  SV2_Results <- SV_summary$BV2_Result
  SV1_Winner <- SV_summary$SV1_Winner
  SV2_Winner <- SV_summary$SV2_Winner
  
  return(c(NSV_Results,Eff_Results,SV1_Results,SV2_Results,SV1_Winner,SV2_Winner))
}


SV_WelfareSummary <- function(tempBootstrap, majorityResults)
{
  # We want to calculate the aggregate welfare of the SV mechanism vs the simple majority
  # Thus, we want to see what the total welfare over the four elections is if we allow
  # the majority to be as in simple majority, versus as if we were in the SV world.
  
  # SV_summary contains the side of the majority with and without BV. Rows 1,2,3,4 correspond
  # to BilingualEduc, Immigration, PublicBonds, and TeachersTenure, respectively. 
  
  # #FOR TESTING ONLY
  # tempBootstrap <- expData
  
  # Aggregate Welfare under Simple Majority
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == SV_summary$NoBV_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == SV_summary$NoBV_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == SV_summary$NoBV_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == SV_summary$NoBV_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_noBV <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  # Aggregate Welfare under Efficiency
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == SV_summary$Eff_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == SV_summary$Eff_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == SV_summary$Eff_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == SV_summary$Eff_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_Eff <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  
  # Aggregate Welfare under Storable Votes - SV1
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == SV_summary$BV1_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == SV_summary$BV1_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == SV_summary$BV1_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == SV_summary$BV1_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_SV1 <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  # Aggregate Welfare under Storable Votes - SV2
  tempRef1 <- filter(tempBootstrap, Ref1_BilingualEduc == SV_summary$BV2_Result[1])
  tempRef1 <- sum(tempRef1$EducationPref)
  tempRef2 <- filter(tempBootstrap, Ref2_TeacherTenure == SV_summary$BV2_Result[4])
  tempRef2 <- sum(tempRef2$TeachersPref)
  tempRef3 <- filter(tempBootstrap, Ref3_PublicBonds == SV_summary$BV2_Result[3])
  tempRef3 <- sum(tempRef3$BondsPref)
  tempRef4 <- filter(tempBootstrap, Ref4_Immigration == SV_summary$BV2_Result[2])
  tempRef4 <- sum(tempRef4$ImmigrationPref)
  
  totalWelfare_SV2 <- sum(tempRef1, tempRef2, tempRef3, tempRef4)
  
  return(c(totalWelfare_noBV, totalWelfare_Eff, totalWelfare_SV1, totalWelfare_SV2))
}

#--- Function to store Vote and Welfare Summaries

saveSummary <- function(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
{
  # Saves what side won in simple majority
  WelfareResults$NSV_Result_BilingualEduc[bootstrapNum] <- voteSummary[1]
  WelfareResults$NSV_Result_Immigration[bootstrapNum] <- voteSummary[2]
  WelfareResults$NSV_Result_PublicBonds[bootstrapNum] <- voteSummary[3]
  WelfareResults$NSV_Result_TeacherTenure[bootstrapNum] <- voteSummary[4]
  
  # # Saves what side should have won under efficient outcomes
  WelfareResults$Eff_Result_BilingualEduc[bootstrapNum] <- voteSummary[5]
  WelfareResults$Eff_Result_Immigration[bootstrapNum] <- voteSummary[6]
  WelfareResults$Eff_Result_PublicBonds[bootstrapNum] <- voteSummary[7]
  WelfareResults$Eff_Result_TeacherTenure[bootstrapNum] <- voteSummary[8]
  
  # Save what side won under SV
  WelfareResults$SV1_Result_BilingualEduc[bootstrapNum] <- voteSummary[9]
  WelfareResults$SV1_Result_Immigration[bootstrapNum] <- voteSummary[10]
  WelfareResults$SV1_Result_PublicBonds[bootstrapNum] <- voteSummary[11]
  WelfareResults$SV1_Result_TeacherTenure[bootstrapNum] <- voteSummary[12]
  
  WelfareResults$SV2_Result_BilingualEduc[bootstrapNum] <- voteSummary[13]
  WelfareResults$SV2_Result_Immigration[bootstrapNum] <- voteSummary[14]
  WelfareResults$SV2_Result_PublicBonds[bootstrapNum] <- voteSummary[15]
  WelfareResults$SV2_Result_TeacherTenure[bootstrapNum] <- voteSummary[16]
  
  # Save whether the minority or majority won (assumes that majority is simple majority)
  WelfareResults$SV1_Winner_BilingualEduc[bootstrapNum] <- voteSummary[17]
  WelfareResults$SV1_Winner_Immigration[bootstrapNum] <- voteSummary[18]
  WelfareResults$SV1_Winner_PublicBonds[bootstrapNum] <- voteSummary[19]
  WelfareResults$SV1_Winner_TeacherTenure[bootstrapNum] <- voteSummary[20]
  
  WelfareResults$SV2_Winner_BilingualEduc[bootstrapNum] <- voteSummary[21]
  WelfareResults$SV2_Winner_Immigration[bootstrapNum] <- voteSummary[22]
  WelfareResults$SV2_Winner_PublicBonds[bootstrapNum] <- voteSummary[23]
  WelfareResults$SV2_Winner_TeacherTenure[bootstrapNum] <- voteSummary[24]
  
  # save welfare summary 
  
  WelfareResults$AggregateWelfare_noBV[bootstrapNum] <- welfareSummary[1]
  WelfareResults$AggregateWelfare_Eff[bootstrapNum] <- welfareSummary[2]
  WelfareResults$AggregateWelfare_SV1[bootstrapNum] <- welfareSummary[3]
  WelfareResults$AggregateWelfare_SV2[bootstrapNum] <- welfareSummary[4]
  
  return(WelfareResults)
}

################################
#### Vote Casting Functions ####
###############################%

#--- casting functions

castPMax <- function(expData, row) # assign vote to a proposal with max points
{
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(expData$EducationPref[row], expData$TeachersPref[row], expData$BondsPref[row], expData$ImmigrationPref[row])
  
  # find index of next max and choose it. 
  nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
  nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
  
  return(viableReferenda[nextRefIndex])
}


castPRand <-  function() # assign vote randomly
{
  randVote <- sample(c("Immigration", "BilingualEduc", "TeacherTenure", "PublicBonds"), 1, prob = rep(.25,4))
  return(randVote)
}


castPEduc <- function(expData, row) # assign vote to an educ proposal with max points
{
  if(expData$EducationPref[row] == expData$TeachersPref[row]) # if both educ prop are equal, randomize
  {
    randVote <- sample(c("BilingualEduc", "TeacherTenure"), 1, prob = c(.5,.5), replace = T)
  } else { # else cast on the one with max points
    maxPointsEduc <- max(expData$EducationPref[row], expData$TeachersPref[row])
    randVote <- ifelse(expData$EducationPref[row] == maxPointsEduc, "BilingualEduc", "TeacherTenure")
  }
  
  return(randVote)
}

#--- check abstain functions

checkVoteOnAbstain_StatModel <- function(expData, row)
{
  #--- first check for SV1
  voteOnAbstain <- T
  
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV1_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV1_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV1_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV1_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      expData$SV1_Strat[row] <- sample(c("pMax", "pRand", "pEduc", "pSal"), 1, prob = SV1_prob, replace = TRUE)
      expData$SV1_Vote[row] <- ifelse(expData$SV1_Strat[row] == "pMax", castPMax(expData, row), 
                                      ifelse(expData$SV1_Strat[row] == "pEduc", castPEduc(expData, row),
                                             ifelse(expData$SV1_Strat[row] == "pRand", castPRand(), "Immigration")))
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  #--- then check for SV2
  voteOnAbstain <- T
  
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV2_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV2_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV2_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV2_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      expData$SV2_Strat[row] <- sample(c("pMax", "pRand", "pEduc", "pSal"), 1, prob = SV2_prob, replace = TRUE)
      expData$SV2_Vote[row] <- ifelse(expData$SV2_Strat[row] == "pMax", castPMax(expData, row), 
                                      ifelse(expData$SV2_Strat[row] == "pEduc", castPEduc(expData, row),
                                             ifelse(expData$SV2_Strat[row] == "pRand", castPRand(), "Immigration")))
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  return(expData[row,])
  
}


checkVoteOnAbstain_Max <- function(expData, row)
{
  #--- first check for SV1
  voteOnAbstain <- T
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(expData$EducationPref[row], expData$TeachersPref[row], expData$BondsPref[row], expData$ImmigrationPref[row])
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV1_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV1_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV1_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV1_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      # find the current issue and eliminate it from consideration
      abstainedRefIndex <- !viableReferenda %in% expData$SV1_Vote[row]
      viableReferenda <- viableReferenda[abstainedRefIndex]
      viablePrefStrength <- viablePrefStrength[abstainedRefIndex]
      
      # find index of next max and choose it. 
      nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
      nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
      
      
      expData$SV1_Vote[row] <- viableReferenda[nextRefIndex]
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  #--- then check for SV2
  voteOnAbstain <- T
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(expData$EducationPref[row], expData$TeachersPref[row], expData$BondsPref[row], expData$ImmigrationPref[row])
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV2_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV2_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV2_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV2_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      # find the current issue and eliminate it from consideration
      abstainedRefIndex <- !viableReferenda %in% expData$SV2_Vote[row]
      viableReferenda <- viableReferenda[abstainedRefIndex]
      viablePrefStrength <- viablePrefStrength[abstainedRefIndex]
      
      # find index of next max and choose it. 
      nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
      nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
      
      
      expData$SV2_Vote[row] <- viableReferenda[nextRefIndex]
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  return(expData[row,])
}


checkVoteOnAbstain_Rand <- function(expData, row)
{
  #--- first check for SV1
  voteOnAbstain <- T
  
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV1_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV1_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV1_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV1_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      expData$SV1_Vote[row] <- castPRand()
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  #--- then check for SV2
  voteOnAbstain <- T
  
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV2_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV2_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV2_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV2_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      expData$SV2_Vote[row] <- castPRand()
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  return(expData[row,])
  
}  


checkVoteOnAbstain_MaxRand <- function(expData, row)
{
  #--- first check for SV1
  voteOnAbstain <- T
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(expData$EducationPref[row], expData$TeachersPref[row], expData$BondsPref[row], expData$ImmigrationPref[row])
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV1_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV1_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV1_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV1_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      
      if(expData$SV1_Strat[row] == "pMax")
      {
        # find the current issue and eliminate it from consideration
        abstainedRefIndex <- !viableReferenda %in% expData$SV1_Vote[row]
        viableReferenda <- viableReferenda[abstainedRefIndex]
        viablePrefStrength <- viablePrefStrength[abstainedRefIndex]
        
        # find index of next max and choose it. 
        nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
        nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
        
        expData$SV1_Vote[row] <- viableReferenda[nextRefIndex]
        
      } else (expData$SV1_Vote[row] <- castPRand())
      
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  
  #--- then check for SV2
  voteOnAbstain <- T
  viableReferenda <- c("BilingualEduc", "TeacherTenure", "PublicBonds", "Immigration")
  viablePrefStrength <- c(expData$EducationPref[row], expData$TeachersPref[row], expData$BondsPref[row], expData$ImmigrationPref[row])
  # check if bv on abstain; if so then sample a new rule and recast vote
  repeat
  {
    if((expData$SV2_Vote[row] == "BilingualEduc" & expData$Ref1_BilingualEduc[row] == "Abstain") |
       (expData$SV2_Vote[row] == "TeacherTenure" & expData$Ref2_TeacherTenure[row] == "Abstain") |
       (expData$SV2_Vote[row] == "PublicBonds" & expData$Ref3_PublicBonds[row] == "Abstain") |
       (expData$SV2_Vote[row] == "Immigration" & expData$Ref4_Immigration[row] == "Abstain"))
    {
      
      if(expData$SV2_Strat[row] == "pMax")
      {
        # find the current issue and eliminate it from consideration
        abstainedRefIndex <- !viableReferenda %in% expData$SV2_Vote[row]
        viableReferenda <- viableReferenda[abstainedRefIndex]
        viablePrefStrength <- viablePrefStrength[abstainedRefIndex]
        
        # find index of next max and choose it. 
        nextRefIndex <- which(viablePrefStrength == max(viablePrefStrength))
        nextRefIndex <- ifelse(length(nextRefIndex) == 1, nextRefIndex, sample(nextRefIndex, 1, prob = rep(1/length(nextRefIndex),length(nextRefIndex))))
        
        expData$SV2_Vote[row] <- viableReferenda[nextRefIndex]
        
      } else (expData$SV2_Vote[row] <- castPRand())
      
    } else {voteOnAbstain <- F}
    
    # if vote has been fixed, then leave loop
    if(voteOnAbstain == F) { break }
  }
  
  return(expData[row,])
  
}  


###################################
##### SV Rule A - AS IS Model #####
##################################%

ruleA <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses the bonus votes as they were cast by the original subjects. 
  # Thus no modification is needed and we just run the welfare measures.
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- 1
  
  voteSummary <- SV_VoteSummary(expData)
  welfareSummary <- SV_WelfareSummary(expData, SV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("SV_RuleAResults", WelfareResults, envir=.GlobalEnv)
}


##################################
##### SV Rule B - Stat Model #####
#################################%

# Reads in the parameters of the statistical model
SV_StatModelData <- read.csv("SV & QV Simulations (No Recalibration)/SV - Raw Simulation Data/SV_StratProb_ParamEstimateMLE.csv", stringsAsFactors = F, header = T)

SV_StatModelData <- summarise(SV_StatModelData, SV1_P_Max = mean(SV1_P_Max), SV1_P_Rand = mean(SV1_P_Rand),
                              SV1_P_Educ = mean(SV1_P_Educ), SV1_P_Sal = mean(SV1_P_Sal), SV2_P_Max = mean(SV2_P_Max), 
                              SV2_P_Rand = mean(SV2_P_Rand), SV2_P_Educ = mean(SV2_P_Educ), SV2_P_Sal = mean(SV2_P_Sal)) 

SV_StatModelData <- round(SV_StatModelData, digits = 2)

SV1_prob <- t(SV_StatModelData)[1:4]
SV2_prob <- t(SV_StatModelData)[5:8]

ruleB <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses the bonus votes as if they occurred using the stat model
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- 1
  
  # Determine the individual subject strategies of casting the BV
  
  expData$SV1_Strat <- sample(c("pMax", "pRand", "pEduc", "pSal"), sampleSize, prob = SV1_prob, replace = TRUE)
  expData$SV2_Strat <- sample(c("pMax", "pRand", "pEduc", "pSal"), sampleSize, prob = SV2_prob, replace = TRUE)
  
  # Implement the strategy
  
  for(row in 1:nrow(expData))
  {
    #### SV1 - Vote cast
    
    expData$SV1_Vote[row] <- ifelse(expData$SV1_Strat[row] == "pMax", castPMax(expData, row), 
                                    ifelse(expData$SV1_Strat[row] == "pEduc", castPEduc(expData, row),
                                           ifelse(expData$SV1_Strat[row] == "pRand", castPRand(), "Immigration")))
    
    #### SV2 - Vote cast 
    
    expData$SV2_Vote[row] <- ifelse(expData$SV2_Strat[row] == "pMax", castPMax(expData, row), 
                                    ifelse(expData$SV2_Strat[row] == "pEduc", castPEduc(expData, row),
                                           ifelse(expData$SV2_Strat[row] == "pRand", castPRand(), "Immigration")))
    
    # check there are no votes on abstain
    expData[row,] <- checkVoteOnAbstain_StatModel(expData, row)
    
  }
  
  voteSummary <- SV_VoteSummary(expData)
  welfareSummary <- SV_WelfareSummary(expData, SV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  assign("SV_RuleBResults", WelfareResults, envir=.GlobalEnv)
}


####################################
##### SV Rule C - BV Max Model #####
###################################%

ruleC <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses the bonus votes on Max
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- 1
  
  
  # Implement the strategy
  
  for(row in 1:nrow(expData))
  {
    #### SV1/SV2 BV on Max Assignment
    
    expData$SV1_Vote[row] <- castPMax(expData,row)
    expData$SV2_Vote[row] <- castPMax(expData,row)
    
    # check no vote is on abstains
    
    expData[row,] <- checkVoteOnAbstain_Max(expData,row)
    
  }  
  
  voteSummary <- SV_VoteSummary(expData)
  welfareSummary <- SV_WelfareSummary(expData, SV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  
  assign("SV_RuleCResults", WelfareResults, envir=.GlobalEnv)
}


######################################
#### SV Rule D - Rand & Max Model ####
#####################################%

ruleD <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses the bonus votes as 50% Random and 50% on Max
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- 1
  
  
  # Determine the individual subject strategies of casting the BV
  
  expData$SV1_Strat <- sample(c("pMax", "pRand"), sampleSize, prob = c(.5,.5), replace = TRUE)
  expData$SV2_Strat <- sample(c("pMax", "pRand"), sampleSize, prob = c(.5,.5), replace = TRUE)
  
  # Implement the strategy
  
  for(row in 1:nrow(expData))
  {
    #### SV1/SV2 BV Assignment 
    
    expData$SV1_Vote[row] <- ifelse(expData$SV1_Strat[row] == "pMax", castPMax(expData,row), castPRand())
    expData$SV2_Vote[row] <- ifelse(expData$SV2_Strat[row] == "pMax", castPMax(expData,row), castPRand())
    
    # check no bv on abstain
    
    expData[row, ] <- checkVoteOnAbstain_MaxRand(expData,row)
  }
  
  
  voteSummary <- SV_VoteSummary(expData)
  welfareSummary <- SV_WelfareSummary(expData, SV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  
  assign("SV_RuleDResults", WelfareResults, envir=.GlobalEnv)
}


######################################
#### SV Rule E - All Random Model ####
#####################################%

ruleE <- function(expData, WelfareResults, bootstrapNum)
{
  # This rule uses the bonus votes as 100% Random
  
  # #FOR TESTING ONLY
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- 1
  
  # Implement the strategy
  
  for(row in 1:nrow(expData))
  {
    #### SV1/SV2 BV Assignment 
    
    expData$SV1_Vote[row] <- castPRand()
    expData$SV2_Vote[row] <- castPRand()
    
    expData[row,] <- checkVoteOnAbstain_Rand(expData,row)
  }
  
  voteSummary <- SV_VoteSummary(expData)
  welfareSummary <- SV_WelfareSummary(expData, SV_summary)
  WelfareResults <- saveSummary(voteSummary, welfareSummary, WelfareResults, bootstrapNum)
  
  
  assign("SV_RuleEResults", WelfareResults, envir=.GlobalEnv)
}


###################################
###### SV Bootstrap Plot Fns ######
##################################%

# summarizes the preferences for a given proposal and puts it in a hist friendly format
# also calculates other preference metrics, e.g., num and total points in favor or against

bootPrefHist <- function(bootstrapData, Proposal)
{
  #---  Creates list of bins for histogram of preferences
  bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                    "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                    "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")
  
  bins10 <- data.frame(Bins = bins_width10, BinOrder = seq(1,21,1), stringsAsFactors = F)
  
  #--- Determines the current proposal
  
  if(Proposal == "BilingualEduc")
  {
    propTitle <- "Bilingual Education"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["EducationPref"]]
  }
  if(Proposal == "Immigration")
  {
    propTitle <- "Immigration"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["ImmigrationPref"]]
  }
  if(Proposal == "TeacherTenure")
  {
    propTitle <- "Teacher Tenure"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["TeachersPref"]]
  }
  if(Proposal == "PublicBonds")
  {
    propTitle <- "Public Bonds"
    bootstrapData[["CurrentPref"]] <- bootstrapData[["BondsPref"]]
  }
  
  ProposalPos <- filter(bootstrapData, CurrentPref > 0) %>% mutate(CurrentPref_Bins = cut(CurrentPref, right = T, breaks = seq(0, 100, 10)))
  ProposalNeg <- filter(bootstrapData, CurrentPref < 0) %>% mutate(CurrentPref_Bins = cut(CurrentPref, right = F, breaks = seq(-100, 0, 10)))
  ProposalZero <- filter(bootstrapData, CurrentPref == 0) %>% mutate(CurrentPref_Bins = "0")
  
  Proposal <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, CurrentPref, CurrentPref_Bins)
  Proposal <- group_by(Proposal, CurrentPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = CurrentPref_Bins)
  Proposal <- full_join(Proposal, bins10, by = "Bins")
  Proposal$Bins <- factor(Proposal$Bins, levels = Proposal$Bins[order(Proposal$BinOrder)])
  
  Proposal[is.na(Proposal)] <- 0
  Proposal <- arrange(Proposal, BinOrder) 
  
  results <- list()
  results[[1]] <- as.data.frame(t(Proposal$NumSubjects))
  results[[2]] <- sum(ProposalPos$CurrentPref)
  results[[3]] <- nrow(ProposalPos)
  results[[4]] <- abs(sum(ProposalNeg$CurrentPref))
  results[[5]] <- nrow(ProposalNeg)
  
  return(results)
}



#############################################
##### Retrieve SV Imm Pref Data #####
############################################%

#--- Get distribution of Imm prefs from population 

# setup the binning dataframe 

bins_width10 <- c("[-100,-90)", "[-90,-80)", "[-80,-70)", "[-70,-60)", "[-60,-50)", "[-50,-40)", 
                  "[-40,-30)", "[-30,-20)", "[-20,-10)", "[-10,0)", "0", "(0,10]", "(10,20]", 
                  "(20,30]", "(30,40]", "(40,50]", "(50,60]", "(60,70]", "(70,80]", "(80,90]", "(90,100]")

bins10 <- data.frame(Bins = bins_width10, BinOrder = seq(1,21,1), stringsAsFactors = F)
bins10 <- mutate(bins10, BinColor = ifelse(Bins=="0",1,0))
binWidth <- 10

# bin and combine population data 

Pop_data <- bind_rows(select(SV_data, MID, matches("Ref")), select(QV_data, MID, matches("Ref")))
# Pop_data <- select(SV_data, MID, matches("Ref")) #***#

ProposalPos <- filter(Pop_data, ImmigrationPref > 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = T, breaks = seq(0, 100, binWidth)))
ProposalNeg <- filter(Pop_data, ImmigrationPref < 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = F, breaks = seq(-100, 0, binWidth)))
ProposalZero <- filter(Pop_data, ImmigrationPref == 0) %>% mutate(ImmigrationPref_Bins = "0")

CurrentProp <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref, ImmigrationPref_Bins)
CurrentProp <- group_by(CurrentProp, ImmigrationPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = ImmigrationPref_Bins)
CurrentProp <- full_join(CurrentProp, bins10, by = "Bins") %>% arrange(BinOrder)
CurrentProp$Bins <- factor(CurrentProp$Bins, levels = CurrentProp$Bins[order(CurrentProp$BinOrder)])
  
# Create vectors with Pop immigration data
ImmPref_Pop <- CurrentProp
ImmPref_Pop[is.na(ImmPref_Pop)] <- 0
binList <- as.character(ImmPref_Pop$Bins) #***#
probDistList <- ImmPref_Pop$NumSubjects #***#

#--- Get distribution of Imm prefs from SV sample 

ProposalPos <- filter(SV_data, ImmigrationPref > 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = T, breaks = seq(0, 100, binWidth)))
ProposalNeg <- filter(SV_data, ImmigrationPref < 0) %>% mutate(ImmigrationPref_Bins = cut(ImmigrationPref, right = F, breaks = seq(-100, 0, binWidth)))
ProposalZero <- filter(SV_data, ImmigrationPref == 0) %>% mutate(ImmigrationPref_Bins = "0")

CurrentProp <- bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref, ImmigrationPref_Bins)
CurrentProp <- group_by(CurrentProp, ImmigrationPref_Bins) %>% summarise(NumSubjects = n()) %>% rename(Bins = ImmigrationPref_Bins)
CurrentProp <- full_join(CurrentProp, bins10, by = "Bins") %>% arrange(BinOrder)
CurrentProp$Bins <- factor(CurrentProp$Bins, levels = CurrentProp$Bins[order(CurrentProp$BinOrder)])

# Create vectors with SV immigration data and MIDs
ImmPref_SV <- CurrentProp
ImmPref_SV[is.na(ImmPref_SV)] <- 0
SV_ImmIDs <-  bind_rows(ProposalNeg, ProposalZero, ProposalPos) %>% select(MID, ImmigrationPref_Bins)


###################################
###### SV Bootstrapping Loop ######
##################################%

bootstrapIterations <- 10000
set.seed(20)
sampleSize <- nrow(SV_data)

# Creates dataframes to store the Rule A-E results

SV_RuleAResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NSV_Result_BilingualEduc = NA, NSV_Result_Immigration = NA, NSV_Result_PublicBonds = NA, NSV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              SV1_Result_BilingualEduc = NA, SV1_Result_Immigration = NA, SV1_Result_PublicBonds = NA, SV1_Result_TeacherTenure = NA, 
                              SV2_Result_BilingualEduc = NA, SV2_Result_Immigration = NA, SV2_Result_PublicBonds = NA, SV2_Result_TeacherTenure = NA,
                              SV1_Winner_BilingualEduc = NA, SV1_Winner_Immigration = NA, SV1_Winner_PublicBonds = NA, SV1_Winner_TeacherTenure = NA, 
                              SV2_Winner_BilingualEduc = NA, SV2_Winner_Immigration = NA, SV2_Winner_PublicBonds = NA, SV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noBV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_SV1 = NA, AggregateWelfare_SV2 = NA)

SV_RuleBResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NSV_Result_BilingualEduc = NA, NSV_Result_Immigration = NA, NSV_Result_PublicBonds = NA, NSV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              SV1_Result_BilingualEduc = NA, SV1_Result_Immigration = NA, SV1_Result_PublicBonds = NA, SV1_Result_TeacherTenure = NA, 
                              SV2_Result_BilingualEduc = NA, SV2_Result_Immigration = NA, SV2_Result_PublicBonds = NA, SV2_Result_TeacherTenure = NA,
                              SV1_Winner_BilingualEduc = NA, SV1_Winner_Immigration = NA, SV1_Winner_PublicBonds = NA, SV1_Winner_TeacherTenure = NA, 
                              SV2_Winner_BilingualEduc = NA, SV2_Winner_Immigration = NA, SV2_Winner_PublicBonds = NA, SV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noBV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_SV1 = NA, AggregateWelfare_SV2 = NA)

SV_RuleCResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NSV_Result_BilingualEduc = NA, NSV_Result_Immigration = NA, NSV_Result_PublicBonds = NA, NSV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              SV1_Result_BilingualEduc = NA, SV1_Result_Immigration = NA, SV1_Result_PublicBonds = NA, SV1_Result_TeacherTenure = NA, 
                              SV2_Result_BilingualEduc = NA, SV2_Result_Immigration = NA, SV2_Result_PublicBonds = NA, SV2_Result_TeacherTenure = NA,
                              SV1_Winner_BilingualEduc = NA, SV1_Winner_Immigration = NA, SV1_Winner_PublicBonds = NA, SV1_Winner_TeacherTenure = NA, 
                              SV2_Winner_BilingualEduc = NA, SV2_Winner_Immigration = NA, SV2_Winner_PublicBonds = NA, SV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noBV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_SV1 = NA, AggregateWelfare_SV2 = NA)

SV_RuleDResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NSV_Result_BilingualEduc = NA, NSV_Result_Immigration = NA, NSV_Result_PublicBonds = NA, NSV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              SV1_Result_BilingualEduc = NA, SV1_Result_Immigration = NA, SV1_Result_PublicBonds = NA, SV1_Result_TeacherTenure = NA, 
                              SV2_Result_BilingualEduc = NA, SV2_Result_Immigration = NA, SV2_Result_PublicBonds = NA, SV2_Result_TeacherTenure = NA,
                              SV1_Winner_BilingualEduc = NA, SV1_Winner_Immigration = NA, SV1_Winner_PublicBonds = NA, SV1_Winner_TeacherTenure = NA, 
                              SV2_Winner_BilingualEduc = NA, SV2_Winner_Immigration = NA, SV2_Winner_PublicBonds = NA, SV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noBV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_SV1 = NA, AggregateWelfare_SV2 = NA)

SV_RuleEResults <- data.frame(Round = seq(1, bootstrapIterations, 1), 
                              NSV_Result_BilingualEduc = NA, NSV_Result_Immigration = NA, NSV_Result_PublicBonds = NA, NSV_Result_TeacherTenure = NA,
                              Eff_Result_BilingualEduc = NA, Eff_Result_Immigration = NA, Eff_Result_PublicBonds = NA, Eff_Result_TeacherTenure = NA,
                              SV1_Result_BilingualEduc = NA, SV1_Result_Immigration = NA, SV1_Result_PublicBonds = NA, SV1_Result_TeacherTenure = NA, 
                              SV2_Result_BilingualEduc = NA, SV2_Result_Immigration = NA, SV2_Result_PublicBonds = NA, SV2_Result_TeacherTenure = NA,
                              SV1_Winner_BilingualEduc = NA, SV1_Winner_Immigration = NA, SV1_Winner_PublicBonds = NA, SV1_Winner_TeacherTenure = NA, 
                              SV2_Winner_BilingualEduc = NA, SV2_Winner_Immigration = NA, SV2_Winner_PublicBonds = NA, SV2_Winner_TeacherTenure = NA,
                              AggregateWelfare_noBV = NA, AggregateWelfare_Eff = NA, AggregateWelfare_SV1 = NA, AggregateWelfare_SV2 = NA)


# Creates dataframes to store ALL preferences for all proposals and bootstraps.
SV_EducBSPref <- data.frame()
SV_BondBSPref <- data.frame()
SV_TeachBSPref <- data.frame()
SV_ImmBSPref <- data.frame()

SV_PrefSummary <- data.frame(Round = seq(1, bootstrapIterations, 1),
                             BilingualEduc_PointsInFavor = NA, BilingualEduc_SubjectsInFavor = NA, 
                             BilingualEduc_PointsAgainst = NA, BilingualEduc_SubjectsAgainst = NA, 
                             Immigration_PointsInFavor = NA, Immigration_SubjectsInFavor = NA, 
                             Immigration_PointsAgainst = NA, Immigration_SubjectsAgainst = NA,
                             PublicBonds_PointsInFavor = NA, PublicBonds_SubjectsInFavor = NA,
                             PublicBonds_PointsAgainst = NA, PublicBonds_SubjectsAgainst = NA,
                             TeacherTenure_PointsInFavor = NA, TeacherTenure_SubjectsInFavor = NA, 
                             TeacherTenure_PointsAgainst = NA, TeacherTenure_SubjectsAgainst = NA)

# Creates dataframes to store the realized individual utilities in each bootstrap
SV_IndivUtility_NSV_A <- data.frame()
SV_IndivUtility_SV1_A <- data.frame()
SV_IndivUtility_SV2_A <- data.frame()

SV_IndivUtility_NSV_B <- data.frame()
SV_IndivUtility_SV1_B <- data.frame()
SV_IndivUtility_SV2_B <- data.frame()

SV_IndivUtility_NSV_C <- data.frame()
SV_IndivUtility_SV1_C <- data.frame()
SV_IndivUtility_SV2_C <- data.frame()

SV_IndivUtility_NSV_D <- data.frame()
SV_IndivUtility_SV1_D <- data.frame()
SV_IndivUtility_SV2_D <- data.frame()

SV_IndivUtility_NSV_E <- data.frame()
SV_IndivUtility_SV1_E <- data.frame()
SV_IndivUtility_SV2_E <- data.frame()

# function to datermine indiv util vectors under diff voting schemes

calcRealizedUtil <- function(BS_data, SV_summary)
{
  SV_temp <- BS_data %>% select(-matches("Abstain")) %>% select(matches("Ref")) %>%
    mutate(Util_NSV_BilingualEduc = ifelse(Ref1_BilingualEduc == SV_summary$NoBV_Result[SV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_SV1_BilingualEduc = ifelse(Ref1_BilingualEduc == SV_summary$BV1_Result[SV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_SV2_BilingualEduc = ifelse(Ref1_BilingualEduc == SV_summary$BV2_Result[SV_summary$Referenda == "BilingualEduc"],EducationPref,0),
           Util_NSV_TeacherTenure = ifelse(Ref2_TeacherTenure == SV_summary$NoBV_Result[SV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_SV1_TeacherTenure = ifelse(Ref2_TeacherTenure == SV_summary$BV1_Result[SV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_SV2_TeacherTenure = ifelse(Ref2_TeacherTenure == SV_summary$BV2_Result[SV_summary$Referenda == "TeacherTenure"],TeachersPref,0),
           Util_NSV_PublicBonds = ifelse(Ref3_PublicBonds == SV_summary$NoBV_Result[SV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_SV1_PublicBonds = ifelse(Ref3_PublicBonds == SV_summary$BV1_Result[SV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_SV2_PublicBonds = ifelse(Ref3_PublicBonds == SV_summary$BV2_Result[SV_summary$Referenda == "PublicBonds"],BondsPref,0),
           Util_NSV_Immigration = ifelse(Ref4_Immigration == SV_summary$NoBV_Result[SV_summary$Referenda == "Immigration"],ImmigrationPref,0),
           Util_SV1_Immigration = ifelse(Ref4_Immigration == SV_summary$BV1_Result[SV_summary$Referenda == "Immigration"],ImmigrationPref,0),
           Util_SV2_Immigration = ifelse(Ref4_Immigration == SV_summary$BV2_Result[SV_summary$Referenda == "Immigration"],ImmigrationPref,0)) %>%
    mutate(Util_NSV = Util_NSV_BilingualEduc + Util_NSV_TeacherTenure + Util_NSV_PublicBonds + Util_NSV_Immigration,
           Util_SV1 = Util_SV1_BilingualEduc + Util_SV1_TeacherTenure + Util_SV1_PublicBonds + Util_SV1_Immigration, 
           Util_SV2 = Util_SV2_BilingualEduc + Util_SV2_TeacherTenure + Util_SV2_PublicBonds + Util_SV2_Immigration)
  
  temp_NSV <- as.data.frame(t(SV_temp$Util_NSV))
  temp_SV1 <- as.data.frame(t(SV_temp$Util_SV1))
  temp_SV2 <- as.data.frame(t(SV_temp$Util_SV2))
  
  return(list(temp_NSV,temp_SV1,temp_SV2))
}




simulatedSamples <- list() # where each simulated dataframe will be stored


# lists to save summary of each simulation data analysis
list_Summary_Pre <- data.frame()
list_Summary_Mid <- data.frame()
list_Summary_Final <- data.frame()

# This loop bootstraps a new sample for SV of the same size using sampling w/ replacement
for(bsCount in 1:bootstrapIterations) 
{
  print(bsCount)
  
  BS_data <- data.frame()
  
  # Use SV_Imm data to create a random sample and then draw from SV_ImmIDs
  randomBinSample <- sample(binList, nrow(SV_data), replace = T, prob = probDistList)
  randomBinSample <- as.data.frame(table(randomBinSample))
  randomBinSample <- rename(randomBinSample, Bin = randomBinSample) %>% mutate(Bin = as.character(Bin))
  
  # Since SV lacks the "[-100,-90)" bin, this makes up the balance with "[-90,-80)"
  randomBinSample$Bin <- gsub("[-100,-90)", "[-90,-80)", randomBinSample$Bin, fixed = TRUE)
  
  for(row in 1:nrow(randomBinSample))
  {
    # Subset the list of names by bin
    SV_temp <- filter(SV_ImmIDs, ImmigrationPref_Bins == randomBinSample$Bin[row])
    
    # Take a random sample and store only the MIDs
    randRowNum <- sample(nrow(SV_temp), randomBinSample$Freq[row], replace = T)
    SV_temp <- SV_temp[randRowNum,]
    SV_temp <- select(SV_temp, MID)
    
    # Merge random MID sample with final MID list
    BS_data <- bind_rows(BS_data, SV_temp)
  }
  
  # Use final bootstrapped list to create a new SV sample
  selectedRows <- match(BS_data$MID, SV_data$MID)
  BS_data <- SV_data[selectedRows,]
  
  # Number the rows to avoid conflicts
  row.names(BS_data) <- seq(1,nrow(BS_data),1)
  
  #save an unedited copy for the preference summary
  BS_data_temp <- BS_data
  simulatedSamples[[bsCount]] <- BS_data
  
  
  # Mutate the data back to all positive & Run the bootstrapped data under the four rules
  BS_data <- mutate(BS_data, ImmigrationPref = abs(ImmigrationPref), BondsPref = abs(BondsPref),
                    EducationPref = abs(EducationPref), TeachersPref = abs(TeachersPref))
  
  # Testing
  # expData <- BS_data
  # WelfareResults <- SV_RuleAResults
  # bootstrapNum <- bsCount
  
  #---------- Run the bootstrapped data under the rules
  
  temp_Pre_All <- data.frame()
  temp_Mid_All <- data.frame()
  temp_Final_All <- data.frame()
  
  ruleA(expData = BS_data, WelfareResults = SV_RuleAResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "A", Testing_SV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "A", Testing_SV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "A", SV_summary))
  utils <- calcRealizedUtil(BS_data, SV_summary)
  SV_IndivUtility_NSV_A <- bind_rows(SV_IndivUtility_NSV_A, utils[[1]])
  SV_IndivUtility_SV1_A <- bind_rows(SV_IndivUtility_SV1_A, utils[[2]]) 
  SV_IndivUtility_SV2_A <- bind_rows(SV_IndivUtility_SV2_A, utils[[3]])
  
  ruleB(expData = BS_data, WelfareResults = SV_RuleBResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "B", Testing_SV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "B", Testing_SV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "B", SV_summary))
  utils <- calcRealizedUtil(BS_data, SV_summary)
  SV_IndivUtility_NSV_B <- bind_rows(SV_IndivUtility_NSV_B, utils[[1]])
  SV_IndivUtility_SV1_B <- bind_rows(SV_IndivUtility_SV1_B, utils[[2]]) 
  SV_IndivUtility_SV2_B <- bind_rows(SV_IndivUtility_SV2_B, utils[[3]])
  
  ruleC(expData = BS_data, WelfareResults = SV_RuleCResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "C", Testing_SV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "C", Testing_SV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "C", SV_summary))
  utils <- calcRealizedUtil(BS_data, SV_summary)
  SV_IndivUtility_NSV_C <- bind_rows(SV_IndivUtility_NSV_C, utils[[1]])
  SV_IndivUtility_SV1_C <- bind_rows(SV_IndivUtility_SV1_C, utils[[2]]) 
  SV_IndivUtility_SV2_C <- bind_rows(SV_IndivUtility_SV2_C, utils[[3]])
  
  ruleD(expData = BS_data, WelfareResults = SV_RuleDResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "D", Testing_SV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "D", Testing_SV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "D", SV_summary))
  utils <- calcRealizedUtil(BS_data, SV_summary)
  SV_IndivUtility_NSV_D <- bind_rows(SV_IndivUtility_NSV_D, utils[[1]])
  SV_IndivUtility_SV1_D <- bind_rows(SV_IndivUtility_SV1_D, utils[[2]]) 
  SV_IndivUtility_SV2_D <- bind_rows(SV_IndivUtility_SV2_D, utils[[3]])
  
  
  ruleE(expData = BS_data, WelfareResults = SV_RuleEResults, bootstrapNum = bsCount)
  
  temp_Pre_All <- rbind(temp_Pre_All, cbind(Round = bsCount, Rule = "E", Testing_SV_Summary_Pre))
  temp_Mid_All <- rbind(temp_Mid_All, cbind(Round = bsCount, Rule = "E", Testing_SV_Summary_Mid))
  temp_Final_All <- rbind(temp_Final_All, cbind(Round = bsCount, Rule = "E", SV_summary))
  utils <- calcRealizedUtil(BS_data, SV_summary)
  SV_IndivUtility_NSV_E <- bind_rows(SV_IndivUtility_NSV_E, utils[[1]])
  SV_IndivUtility_SV1_E <- bind_rows(SV_IndivUtility_SV1_E, utils[[2]]) 
  SV_IndivUtility_SV2_E <- bind_rows(SV_IndivUtility_SV2_E, utils[[3]])
  
  # combine and save data 
  
  list_Summary_Pre <- rbind(list_Summary_Pre, temp_Pre_All)
  list_Summary_Mid <- rbind(list_Summary_Mid, temp_Mid_All)
  list_Summary_Final <- rbind(list_Summary_Final, temp_Final_All)
  
  

  
  #--- Preference Histogram and other Preference Summaries 
  
  # BilingualEduc
  PrefHistResults <- bootPrefHist(BS_data_temp, "BilingualEduc")
  SV_EducBSPref <- bind_rows(SV_EducBSPref, PrefHistResults[[1]])
  
  SV_PrefSummary$BilingualEduc_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  SV_PrefSummary$BilingualEduc_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  SV_PrefSummary$BilingualEduc_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  SV_PrefSummary$BilingualEduc_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  
  # Immigration 
  PrefHistResults <- bootPrefHist(BS_data_temp, "Immigration")
  SV_ImmBSPref <- bind_rows(SV_ImmBSPref, PrefHistResults[[1]])
  
  SV_PrefSummary$Immigration_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  SV_PrefSummary$Immigration_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  SV_PrefSummary$Immigration_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  SV_PrefSummary$Immigration_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  # TeacherTenure 
  PrefHistResults <- bootPrefHist(BS_data_temp, "TeacherTenure")
  SV_TeachBSPref <- bind_rows(SV_TeachBSPref, PrefHistResults[[1]])
  
  SV_PrefSummary$TeacherTenure_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  SV_PrefSummary$TeacherTenure_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  SV_PrefSummary$TeacherTenure_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  SV_PrefSummary$TeacherTenure_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
  
  # PublicBonds 
  PrefHistResults <- bootPrefHist(BS_data_temp, "PublicBonds")
  SV_BondBSPref <- bind_rows(SV_BondBSPref, PrefHistResults[[1]])
  
  SV_PrefSummary$PublicBonds_PointsInFavor[bsCount] <- PrefHistResults[[2]]
  SV_PrefSummary$PublicBonds_PointsAgainst[bsCount] <- PrefHistResults[[4]]
  
  SV_PrefSummary$PublicBonds_SubjectsInFavor[bsCount] <- PrefHistResults[[3]]
  SV_PrefSummary$PublicBonds_SubjectsAgainst[bsCount] <- PrefHistResults[[5]]
}


# Rename the column names of the dataframes with the preferences for each bootstrap and each referenda
names(SV_EducBSPref) <- gsub("V","EducPref",names(SV_EducBSPref))
names(SV_ImmBSPref) <- gsub("V","ImmPref",names(SV_ImmBSPref))
names(SV_BondBSPref) <- gsub("V","BondPref",names(SV_BondBSPref))
names(SV_TeachBSPref) <- gsub("V","TeachPref",names(SV_TeachBSPref))


# Create a column to store the Round number for the list of individual utilities for each bootstrap
SV_IndivUtility_NSV_A <- tibble::rownames_to_column(SV_IndivUtility_NSV_A, "Round")
SV_IndivUtility_NSV_B <- tibble::rownames_to_column(SV_IndivUtility_NSV_B, "Round")
SV_IndivUtility_NSV_C <- tibble::rownames_to_column(SV_IndivUtility_NSV_C, "Round")
SV_IndivUtility_NSV_D <- tibble::rownames_to_column(SV_IndivUtility_NSV_D, "Round")
SV_IndivUtility_NSV_E <- tibble::rownames_to_column(SV_IndivUtility_NSV_E, "Round")
SV_IndivUtility_SV1_A <- tibble::rownames_to_column(SV_IndivUtility_SV1_A, "Round")
SV_IndivUtility_SV2_A <- tibble::rownames_to_column(SV_IndivUtility_SV2_A, "Round")
SV_IndivUtility_SV1_B <- tibble::rownames_to_column(SV_IndivUtility_SV1_B, "Round")
SV_IndivUtility_SV2_B <- tibble::rownames_to_column(SV_IndivUtility_SV2_B, "Round")
SV_IndivUtility_SV1_C <- tibble::rownames_to_column(SV_IndivUtility_SV1_C, "Round")
SV_IndivUtility_SV2_C <- tibble::rownames_to_column(SV_IndivUtility_SV2_C, "Round")
SV_IndivUtility_SV1_D <- tibble::rownames_to_column(SV_IndivUtility_SV1_D, "Round")
SV_IndivUtility_SV2_D <- tibble::rownames_to_column(SV_IndivUtility_SV2_D, "Round")
SV_IndivUtility_SV1_E <- tibble::rownames_to_column(SV_IndivUtility_SV1_E, "Round")
SV_IndivUtility_SV2_E <- tibble::rownames_to_column(SV_IndivUtility_SV2_E, "Round")


# Save and export the dataframes with the welfare results and bootsrapped preferences 
write.csv(SV_EducBSPref, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_EducBSPref.csv", row.names = F)
write.csv(SV_ImmBSPref, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_ImmBSPref.csv", row.names = F)
write.csv(SV_BondBSPref, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_BondBSPref.csv", row.names = F)
write.csv(SV_TeachBSPref, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_TeachBSPref.csv", row.names = F)

write.csv(SV_PrefSummary, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_PrefPointSummary.csv", row.names = F)
write.csv(SV_IndivUtility_NSV_A, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_A.csv", row.names = F)
write.csv(SV_IndivUtility_NSV_B, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_B.csv", row.names = F)
write.csv(SV_IndivUtility_NSV_C, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_C.csv", row.names = F)
write.csv(SV_IndivUtility_NSV_D, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_D.csv", row.names = F)
write.csv(SV_IndivUtility_NSV_E, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_NSV_E.csv", row.names = F)
write.csv(SV_IndivUtility_SV1_A, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_A.csv", row.names = F)
write.csv(SV_IndivUtility_SV2_A, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_A.csv", row.names = F)
write.csv(SV_IndivUtility_SV1_B, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_B.csv", row.names = F)
write.csv(SV_IndivUtility_SV2_B, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_B.csv", row.names = F)
write.csv(SV_IndivUtility_SV1_C, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_C.csv", row.names = F)
write.csv(SV_IndivUtility_SV2_C, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_C.csv", row.names = F)
write.csv(SV_IndivUtility_SV1_D, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_D.csv", row.names = F)
write.csv(SV_IndivUtility_SV2_D, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_D.csv", row.names = F)
write.csv(SV_IndivUtility_SV1_E, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV1_E.csv", row.names = F)
write.csv(SV_IndivUtility_SV2_E, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_IndivUtility_SV2_E.csv", row.names = F)

write.csv(list_Summary_Pre, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_Summary_Pre.csv", row.names = F)
write.csv(list_Summary_Mid, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_Summary_Mid.csv", row.names = F)
write.csv(list_Summary_Final, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_Summary_Final.csv", row.names = F)


write.csv(SV_RuleAResults, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleAResults.csv", row.names = F)
write.csv(SV_RuleBResults, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleBResults.csv", row.names = F)
write.csv(SV_RuleCResults, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleCResults.csv", row.names = F)
write.csv(SV_RuleDResults, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleDResults.csv", row.names = F)
write.csv(SV_RuleEResults, "./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/SV_RuleEResults.csv", row.names = F)

saveRDS(simulatedSamples, file="./SV & QV Recalibration on Pop Imm/SV - Raw Simulation Data/simulatedSamples.RData")
